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Abstract 

A model for a flat isotropic universe with a negative cosmological constant A and a massless 
scalar field as sole matter content is studied within the framework of Loop Quantum Cosmology. 
By application of the methods introduced for the model with A = 0, the physical Hilbert space and 
the set of Dirac observables are constructed. As in that case, the scalar field plays here the role of 
an emergent time. The properties of the system are found to be similar to those of the k = 1 FRW 
model: for small energy densities, the quantum dynamics reproduces the classical one, whereas, 
due to modifications at near-Planckian densities, the big bang and big crunch singularities are 
replaced by a quantum bounce connecting deterministically the large semiclassical epochs. Thus 
in Loop Quantum Cosmology the evolution is qualitatively cyclic. 
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I. INTRODUCTION 



Loop Quantum Cosmology [1] - an application of methods of Loop Quantum Grav- 
ity jl] to symmetry reduced models - constitutes a promising way of studying quantum- 
gravitational effects in cosmological models. In particular one of the simplest models, a flat 
Friedmann- Robertson- Walker (FRW) universe was analyzed within its framework 0, 0, 0] ■ 
In that case, the structure of the Hamiltonian constraint allowed to treat the constrained 
system as a free one, evolving with respect to the scalar field which thus plays the role of 
an emergent time. This, in turn, allowed the construction of a physical Hilbert space and a 
set of Dirac observables, which were used next to extract the physics by means of numerical 
methods. The results were quite surprising: the analysis has shown that, when the matter 
energy density approaches the Planck scale, the quantum-geometric effects cause gravity to 
become repulsive. In consequence, a large semiclassical expanding universe is preceded by a 
(also large and semiclassical) contracting one, deterministically connected to the former by 
a quantum bridge. The transition point of the evolution (called quantum bounce) is char- 
acterized by an energy density which, at this point, equals the critical value p c ~ 0.82p P1 . 
Furthermore, even when quantum corrections actually dominate the dynamics, the state 
representing the universe remains semiclassical - its evolution is to great precision described 
by the so called classical effective dynamics (1, @] . 

The results obtained for the flat FRW model were next generalized to the spherical one [7] 
(the k = 1 FRW model). The properties of the Hilbert space and an evolution operator were 
investigated analytically (1, 0] and the robustness of their features was confirmed through 



the analysis of its approximation (known as sLQC) [10l . Further generalizations to 

anisotropic (and further inhomogeneous) models by different research groups are in various 
stages of progress [H, 13. lH . 



Thus far, however, the only models described rigorously were universes with a vanishing 
cosmological constant A and a massless scalar field. In this article, we extend the analysis 
of [BJ to include the universes with negative A. Although the observations favor a posi- 
tive A, this model constitutes a convenient way of testing which features of the previously 
investigated model we can hope to generalize to more realistic systems. Also, since it is a 
classically recollapsing system, we can use it to investigate semiclassicality issues (dispersion 
after many 'cycles' of evolution). The specific questions we intend to address here are the 
following: 

• Do the qualitative features of the A = model survive also in this case? In particular, 
are the big bang/crunch singularities replaced by quantum bounces as in the previously 
investigated cases? All the models analyzed so far not only experienced the bounce, 
but for Gaussian states the observed dispersion of the wave packet after the bounce 
was severely restricted by the values of the spreads before it. In the flat case this 
result was next generalized analytically to a space of states admitting semiclassical 
epoch 1 [ll[ within the context of sLQC. Therefore, it is important to ask whether 
such behavior will occur also in the considered model, or it was just a result of the 
extreme simplicity of the previous ones. 

• If the answer to the previous question is in the affirmative, then is the critical energy 



The states for which either at early or late times the relative dispersions of chosen Dirac observables are 
< 1. 
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density p c still a fundamental bound? In both the k = and k = 1 models for 
physically sensible 2 states, the matter energy density at the bounce point agreed to 
great precision with p c . Furthermore, later investigations within the sLQC model have 
shown that p c is indeed a fundamental energy bound. But again, we do not know a 
priori whether this feature is characteristic just of the models investigated so far and 
how (if at all) it generalizes. 

• Does this model possess any new feature not observed in A = or k = 1 case ? 

A preliminary investigation of the A < model has been conducted already in However, 
the physical Hilbert space was not constructed; the goal there was only to verify the per- 
sistence of the bounce. Recently, a heuristically constructed effective classical Hamiltonian 
was used [l5| to obtain the effective trajectories of both the A < and A > systems and 
analyze the effect of the quantum-geometric corrections on the universe's dynamics. How- 
ever, since the effective Hamiltonian was not derived systematically, the results have to be 
confirmed against genuine quantum evolution. 

In addition to the problems described above, we also address the concerns about the choice 
of the symmetric sector of the physical Hilbert space that is sometimes raised. Because of the 
absence of fermions, the triad orientation reflection is a large gauge symmetry. This allowed 
one to restrict the physical Hilbert space to the states symmetric under parity reflection. 
However, since the choice of antisymmetric states is equally justified, it is natural to ask 
whether the results of LQC are robust and will continue to hold if the antisymmetric sector 
is chosen. We address this issue by analyzing, in addition to the standard symmetric states, 
also the space of antisymmetric ones and establish robustness. 

The paper is organized as follows: we start with a brief summary of the basic framework 
(introduced already in earlier papers) in section [III Its content is divided into three parts: 
the classical theory, the kinematics of LQC and the derivation of the quantum Hamiltonian 
constraint. In section II I II we consider a geometrodynamical equivalent of the model - the 
Wheeler-DeWitt (WDW) one. The reason for that is two- fold: first, it will allow us to com- 
pare the results of LQC against a standard quantum model and identify the nonperturbative 
quantum-geometric effects. Second, it will serve as an introduction to the methodology of 
extracting physics, used next on the LQC model. The analytical solvability of the WDW 
model will allow us to show these methods without having to deal with the complications 
of numerical analysis. Analysis of the physical sector is carried out in section IIVI There, 
we extensively use the results of the numerical study described in turn in section [V] That 
section contains also a description of the construction and analysis of the states semiclassical 
at late times. The final results and their discussion are placed in section IVII 

Apart from the main body, the article contains two appendices: in section |XJ we analyze 
the space of antisymmetric states, whereas IBI contains a description of the heuristic methods 
used to extract some of the results. 



2 This indicates the states of the scalar field with momentum sufficiently high for the closed universe to 
grow to macroscopic (> 1 megaparsec) scales before recollapsing. 
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II. THE LQC QUANTIZATION SCHEME 

In this section, we introduce the quantization framework used in later sections of the 
paper. Since we directly apply the framework described in detail in (1, 0], we will just 
present a brief sketch of it. For a more detailed discussion, the reader is referred to the 
above mentioned articles. 

The content of this section is divided into three parts. In the first, we present the classical 
theory used as a basis for quantization. The second part is dedicated to the description of 
the LQC kinematics. Finally, we recall the derivation of the LQC Hamiltonian constraint. 



A. Classical theory 

A flat (k = 0) FRW model represents a spacetime admitting a foliation by spatial isotropic 
3-surfaces M of topology M 3 . Its metric tensor can be written in the form 

g = -dt 2 + a 2 (t)°q , (2.1) 

where t is a time parameter (the cosmic time), °q is a unit (fiducial) Cartesian metric on 
the surface M and the function a(t) is called a scale factor. 

Due to the homogeneity and noncompactness of M, one cannot write an action or Hamil- 
tonian as an integral of the appropriate density over the entire M. Instead, we can define 
them as integrals over a chosen fiducial cubical cell V, constant in comoving coordinates 3 . 
Given such a cell, one can define a triad °e (and cotriad to dual to it) as directed along the 
edges of V and orthonormal with respect to °q. 

As gravitational phase space variables, we choose the connections A\ and the density- 
weighted triads 
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A\ = cV~ « °u l a , E? = pV X . O 

where V Q is a volume of V with respect to °q. The real parameters c, p called respectively 
connection and triad coefficients coordinatize the (2-dimensional) phase space of the gravita- 
tional degrees of freedom. Appropriate scaling by V ensures the invariance of the symplectic 
structure of this phase space (when written in terms of c,p) under different choices of °q. 
The Poisson bracket between c and p equals 

{c,p} = — ^— , (2.3) 

where 7 is the Barbero-Immirzi parameter. 

The basic variables defined as in ( 12.21) automatically satisfy the Gauss and diffeomor- 
phism constraints. The contribution of the geometry to the only nontrivial constraint - the 
Hamiltonian one - is of the form 

C grav = -\ I d 3 x(e ijk e- 1 E ai E b SF5 > ->fA) = -^c^ + Apl (2.4) 
7 *j v 7 



3 The considered model is of the Bianchi type A: the equations of motion derived from the Hamiltonian 
specified in this way are identical to the Einstein field equations reduced to the isotropic case. 
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where e := y/\ detE\ and the field strength F* b := 2d a A b t + e^A l a A j b . 

The only matter content - a homogeneous massless scalar field - is described by two 
global variables: the field value and its conjugate momentum p^, with Poisson bracket 
between them 

{<f>,P4>} = 1 • (2-5) 

The pair (4>,p<p) coordinatizes the phase space corresponding to the matter degrees of free- 
dom. The full phase space of the system is thus 4- dimensional. The complete Hamiltonian 
constraint is of the form 
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C =: C gvSLV + C,p = , where C+ = 9mGp~^ . (2.6) 

The above constraint defines a 3D hypersurface in the 4D phase space. Furthermore, since 
C does not depend explicitly on 0, the momentum p^ is a constant of motion. Therefore, 
the dynamical trajectories can be represented as a (parametrized by p^) family of functions 
p{4>) 

P W = 777 1 W /To — 77/ i TTT ( 2 ' 7 ) 

| A| a cosh(V127rG(0 - O )) 

Their form implies that the considered system recollapses. Each trajectory starts at the big 
bang singularity and ends in a big crunch. 



B. Kinematics of LQC 

To quantize the system, we follow the Dirac program. First we construct a kinematical 
Hilbert space: in our case, it is the tensor product of spaces corresponding to, respectively, 
gravitational and matter degrees of freedom: Tt km = 7^g™ v ® Ti.^ m . 

For the matter we apply the standard Schrodinger quantization. As 'H ] ^ n we choose the 
standard Hilbert space of square integrable functions H^ m = L 2 (IR, d0). The basic operators 
are and p^. To describe the state we choose the (dual) basis (0| of eigenstates of 0. The 
action of 0, p^ on the state can be then expressed as follows 

0^(0) = 0^(0), p^>{4>) = -ihd^icj)) , where ^(0) := (0|^) . flU} 

The quantization of the gravitational degrees of freedom within LQC at the kinematical 
level has been rigorously performed in [16j. The procedure is the analog of the quanti- 
zation scheme used in full LQG (see for example 13] )• Here the basic variables are tri- 
ads and connections along straight edges generated by °ef . The kinematical Hilbert space 
is the space of square integrable functions on the Bohr compact ification of the real line 
^gmv = -^ 2 (^Bohr, d/iBohr)- We will represent its elements using the basis consisting of the 
eigenf unctions of p (promoted to an operator), labeled by 6 1. Despite the continuity of 

the elements of the chosen basis are orthonormal with respect to Kronecker delta 

(/ii|/i 2 ) = <W • (2.9) 
As basic quantum operators, we select p and exp(«4f ) 4 . Their action on the basis elements 



Since the family exp(iAc/2) is not weakly continuous, the operator c does not exists. 



5 



is given by: 

m = , Sq5(^)|/i) = |M + A) ■ dSUD 

Since the holonomy along the edge of fiducial length A generated by °e° can be expressed 
via exp(iAc/2) 

h[ X) = ±[exp(^)+exp(-^)]I+i[exp(^)-exp(-^)]r fc (2.11) 

(where the are related to the Pauli matrices via 2ir k — its quantum analog h!^ 
can be expressed in terms of the operators exp in the same way. 



C. LQC: the Hamiltonian constraint 

In order to write the quantum operator corresponding to the Hamiltonian constraint 
TM, 12.61) . we need to reexpress it in terms of the basic objects selected in the previous 
subsection. 

Let us start with C grav (12.41) . The quantization of the cosmological term is straightforward 
(and just amounts to promoting p to operator p). The remaining part is an integral of the 
product of two terms: e~ x E a% E h i and F^ b . 

Following Thiemann [18j|, we can rewrite the first term in the following form 

e ijk e- l E ai E^ = V Sgn(p) t °e abc °u k c Tr (h™{h£ r \ V}ri) (2.12) 
* 2tt 7 GAK 3 V J 

where V = |pp is the (physical) volume of the cell V. 

The field strength term F^ b can, on the other hand, be approximated via holonomies 
along the square loop oriented on the i-j plane. 

fh [x) - l\ 

= -2 Tr f ) r k ^ ^ , = h<Ph™h<*> . « 



X 2 Vo a 



The size of is fixed by the requirement that its physical area equals the lowest nonzero 
eigenvalue of the LQG area operator 

A = Jl{p) s.t. Ar Di . = f\p\ = A := (2^7)^ . (2.14) 

To express the action of the operator corresponding to it is convenient to use, instead 
of the label /i, a new label v defined as follows 

3 2\/2 

v := Kaga(ji)\fi\2 , K := . (2.15) 

3^3^ 

In the new labeling an exponent operator exp(^) -the component of (via (12. lip )- acts 
simply as a unit translation 

e£p(ific)\v) = \v + l) . (2.16) 
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In the matter part of the Hamiltonian constraint, the only nontrivial component is \p\ 3 / 2 , 
but again this can be reexpressed in terms of holonomies via Thiemann's method 



sgn(» 



Tr^T*hf{hp-\vi} 



[2.17) 



Finally, applying all the results (12.11112.171) to (12.61) . one can write the operator C. We do 
so choosing, in the process, a particular factor ordering (the so called Kaminski ordering) 
0, in which C grav is manifestly symmetric and positive-definite. The action of the final 
result on the state ^> G 7i km can be written in the following form 



where ty(v,<fi) := (v,cj)\^f) and the functions B(v), C\(v) equal 



B(v) :- 



27 K 



\v + 1| 3 — \v — 1| 3 



C A (v) :-- 



167T 2 7 3 4l AI 



and B is an operator corresponding to the A = case derived in [5| 

e o V(v,(P) = -[B(v)]- 1 (C + (v)^(v + 4,^) + C (v)^(v,(P) + C^(v)^(v-4, ( p)) 

with coefficients C7 ± , C°, equal to 
3irKG , 



C + {v) 
C-(v) 



2 v+3 - v 



C + (v-A) 



C°(v) 



-C + {v)-C-{v) . 



(2.18) 

(2.20) 

(2.21a) 
(2.21b) 



For reasons we will explain in later sections of the paper, the operator is called an evolution 
operator. It is symmetric and positive-definite (with respect to the measure B(v )d/iBohr) on 
the domain T> of finite linear combination of states \v). 



III. THE WHEELER-DEWITT LIMIT 



The quantization scheme presented in the previous section is motivated by LQG; however, 
it is not the only method applicable to the system. By replacing Ti^^ with 21g™ v := 
L 2 (R,dfi) and taking the limit A — » in expressions (12 . 1 21 12.131 12 . 1 71) . one arrives to the 
system equivalent to the one originating from geometrodynamics, known as the Wheeler- 
DeWitt system. In the literature, the system obtained from LQC via this procedure is called 
a WDW limit. We will study it in this section in order to identify the effects of the spacetime 
discreteness. We will keep this terminology in the paper although (as it was shown in [10]) 
the WDW model is not the limit of the LQC model in any precise sense. One should think 
about it as the WDW equivalent of an LQC model. 



A. WDW constraint equation, emergent time 



The evolution operator is a sum of two terms: a A = operator O and a A-dependent 
potential term (12.181) . The WDW limit of O was derived in [5( and is of the form 

© *M) = -127rG(^) 2 ^,0) (3.1) 
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where \j/ G 7i km := 21grav ® Calculating the limit of the cosmological constant term 

requires just replacing B in the potential term by its point limit B_ := K/\v\ for A — > 0. In 
consequence, the WDW equivalent of equation (12.1 8p has the form 

dl 0) = -6 0) = 12nG(vd v ) 2 mv, 0) + ijK^h ^-^ ^ ' (3 ' 2) 

where the operator is symmetric and positive-definite with respect to the measure Bdv 
in the standard domain of fast-decaying functions (Schwartz space). 

The above constraint divides the domain of v into two independent sectors, corresponding 
to different signs of v, i.e. to different orientations of the triad Ef. Due to the absence of 
a parity violating interaction in the considered system, we can restrict the studies to states 
that are symmetric/antisymmetric with respect to a reflection in v. For further analysis, we 
choose the symmetric sector, that is \j/(0, v) = —v); however, the presented construction 
can be repeated directly also in the antisymmetric case, with equivalent results. 



B. General solutions, frequency decomposition 

The constraint (I3.2p is similar in its form to the Klein-Gordon equation. Furthermore, 
since there is no explicit dependence on in either (12.61) or (13.21) . is a constant of motion 
of both the classical and the quantum system. Also, at the classical level <fi is monotonic in 
time: we can thus follow the prescription of 0] and reinterpret the constraint, treating it as 
an evolution equation of a free system evolving with respect to <p. The scalar field becomes 
then an emergent time as in the case A = 0. 

To construct the physical Hilbert space we need to find the spectrum of the self-adjoint 
extension of O. The eigenfunction corresponding to an eigenvalue u 2 satisfying 

Ifi-TT 2 ^ 3 / 4 

u 2 ±{v) = -12nG(vd v ) 2 ±(v) - 27Km Av2 ^ v) (3 ' 3) 
can be written in terms of Bessel functions of the third kind 

Mv) = c {I) l ik (f3V^A\v\) + c w /C ifc ( / 3 v /z A>|) , (3.4) 

where k := u/Vl2irG, := 2^/^fh£ P i/{9K) and c {I) ,c (K) G C. When pV^X\v\ < k, 
both I and /C show oscillatory behavior. In particular, as \v\ — > 0, they approach the 
eigenfunctions of the B G operator corresponding to the same frequency uj 

4>uj(v) = c + exp(ik In \v |) + c~ exp(— ik In \v |) . (3.5) 

The complex coefficients c + , c~ of the limit can be determined uniquely as functions of C(j), 

C{K)- 

For j3^/— A\v | > k, the functions X grow exponentially, whereas the functions K expo- 
nentially decay. In consequence, only the eigenfunctions with cm = will contribute to the 
spectral decomposition of 6. This implies that the spectrum of 6 equals Sp(G) = [0, oo) 
and is continuous. Furthermore, due to (13. 5ft . the eigenfunctions with cm = are Dirac 
delta normalizable. Therefore, we can choose the basis setting := a(u)JCik(PV— M v \)i 
where a is a real, positive, a;-dependent normalization factor chosen to satisfy the relation 

(eJs,) = S(u,u') . (3.6) 
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At this point, we note that the structure of the spectral decomposition of is similar 
to the one of the WDW limit for the k — 1 FRW model so that we can follow the 
construction used there. Each element ip{v) of L 2 (M., B_(v)dv ) can be decomposed in the 
basis e^: 

POO 

ij>(v) = / du^(u;)e» . (3.7) 
Jo 

where ip G L 2 (R, du). Therefore, the solutions to the evolution equation (13. 2p with initial 
data in the Schwartz space can be represented in terms of the two functions ^/±(uj) G 
L 2 (R,du) 

*(v,<f>) = Jdu;[* + {u)^(v)e^ + *4u)z a {v)e- iu +\ . (3.8) 

The solutions with vanishing £ + and (denoted in the following as are called the 

negative and positive frequency solutions respectively. Their general form can be written in 
terms of the square root of the operator; namely, for initial data ip (v) specified at (ft = <p Q , 
we have: 

V ± (v,(f>) = e^^-^tpoiv) . (3.9) 



C. Physical Hilbert space, observables 

To construct the physical Hilbert space 7i phy , we again follow (1, 0]. As is the sum 
of the O operator (which is just df n ^) and the positive potential term, it is essentially 
self-adjoint and positive-definite Friedrich's extension of it is thus a unique self-adjoint 
one. One can then apply group averaging techniques [l9| (see the discussion in j3]) to find 
H phy and the inner product. The result is the following: the space 7i phy itself consists of 
normalizable solutions to (j3.2P ; however, as the spaces of positive and negative frequency 
solutions are superselected sectors, we can take as 7i phy the restriction to just one of them. 
Following previous works, we chose positive frequency part, thus defining 7i phy as: 

My A) = J du§L{u)e iJ e i0J ' t> , £ G L 2 (R+ , du) . (3.10) 

The physical inner product within this space can be written as: 



(_vj/|$) = / B(v)dvV(y)$(v) . (3.11) 

J (f>=4> 

In order to be able to extract physical information out of our system, we need to define 
a set of Dirac observables, i.e. self-adjoint operators preserving H phy . Here again we can 
directly use the scalar field momentum and , the amplitude of v at a given 0, defined 
already for A = and k — 1. Their action on the elements \l/ of 7i phy is the following 

p^-qj = -ihd^ , = eV&(*-*'>|i;|tf(v J / ) , < KT3t 

and their expectation values equal respectively: 

&M) = ~ ih f B(v)dv§L(v,(j))(d <p mv,<i>) , (3.13a) 

J </)=const 

= J B(v)dv\v\\^L(v,<p)\ 2 . (3.13b) 
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D. Semiclassical states 



Once we have the physical Hilbert space, the inner product and the observables, we can 
investigate the evolution of a universe represented by a given state. A particularly interesting 
question one can ask is whether, in the considered system, the singularity is resolved. To 
address this question, we construct a Gaussian state which, at a given time <f) , is sharply 
peaked at a large scalar field momentum = hw* (with spread a/y/2) and volume v* and 
is expanding: 



du e 



iv e 



lLj(<j)—<j>*) 



where 



1 



: arcosh 



/ 3Ky/l2irG P% 
{ (4VT740 372 J\A\v 



(3.14) 



(3.15) 



Because of the complicated form of e^, the wave function (13.141) and expectation values 
(I3.13P were calculated numerically (see section |V] for the details). An example of the results 
is shown on Figs [H The state remains semiclassical (sharply peaked) and simply follows the 
classical trajectory (12.71) 



v{cj>) 



(4^7^)1 



cosh (Vi2ttG(0 - 



to the big bang and big crunch singularities, 
the classical singularities are not resolved. 



Wv,«| 



a) 




(3.16) 



In consequence, similarly to the A = case, 



4> 



b) 



WDW 
classical 



10 J 
|v| 



FIG. 1: An example of a Wheeler-DeWitt Gaussian wave packet (|3.14p generated for the parameter 
values A = —0.01, = 5 ■ 10 3 , Ap^/p^ = 0.02 and 4>* = 0. Fig.a) shows the absolute value of the 
wave function. For the presentation clarity, only the points of |\P(u, (f>)\ > 10 -6 were plotted. Fig. b) 
presents the expectation values and dispersions of |£>|<^ (red bars) compared against the classical 
trajectory v(<p) (blue line). As we can see, the quantum trajectory agrees with the classical one 
(the difference being much smaller that the spread). Due to the large changes in magnitude of v 
during the evolution, the trajetory was plotted in logarythmic scale. 
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IV. PHYSICAL SECTOR OF LQC 



The analysis in the previous section allowed to find dynamics predicted by the WDW 
limit of the considered LQC model. Now we perform an analogous study of the model of 
interest. Due to qualitative similarities of the Hamiltonian constraint with its WDW limit, 
the analysis can be performed analogously to the one done in section II II I (with only slight 
modifications required by the fact that is now a difference operator). Following that work, 
we again restrict the study to states symmetric under parity reflection. 5 

First we note that, thanks to the fact that G is a difference operator, we can natu- 
rally divide the gravitational kinematical Hilbert space onto superselected sectors 7^" v = 
©ee[o 2] ^grav,e> where 7^grav,e are the restrictions of 7i.g^ v to the functions supported on the 
sets C e := {±e + An; jigZ} preserved by the action of the Hamiltonian constraint (I2.18P 
and parity reflection LT : ip(v) i— > v). Following the literature, we call these sets lat- 
tices and work with sine le sectors H kin £ := W£w,e ® ?^ in - The kinematical inner product 
corresponding to them is just a restriction of the product of 7i km . 

For each of the sectors illustrated above, the operator 6 is obviously well defined and 
symmetric (with respect to the measure S(v)d/iBohr) on the domain D e - the space of finite 
combinations of \v) with v G C e . Its mathematical properties were rigorously analyzed in 
[§]. It is essentially self-adjoint, its extension is positive-definite and its spectrum is discrete. 
The first two properties allow us again to choose as an emergent time and treat as an 
evolution operator. 

The discreteness of G's spectrum implies that the eigenf unctions relevant for its spectral 
decomposition are normalizable. Furthermore, a numerical study (discussed in section |V|) 
shows that the spectrum is nondegenerate. In consequence, for each allowed value of the 
label s, we can build the physical Hilbert space 7i phj e as a space of normalizable positive 
frequency solutions to (12.181) , analogously to the construction in sections IIIIBI and IIII CI 

zfytf = Ve*, = ^* n e»e^ , (4.1) 

neN 

where ^ are square summable and e n (v) are symmetric in v and normalized eigenfunctions 
of 6, corresponding to eigenvalues which form the basis of 7i phy e . The physical inner 
product can be found through group averaging analogously to the WDW case and can be 
written in the form 

oo 

= J]$A = B{v)*{v)$(v) . (4.2) 

n=0 veC e 

To complete the quantization program we need to choose a set of Dirac observables. 
In order to be able to compare the results with the WDW limit, we choose the operators 
analogous to (I3.12p 

= -%nd^l , = e iv ^-*' ) |u|tf(u,0 / ) ■ (Q 



5 It is also correct to work with the antisymmetric sector of the theory. We discuss that case in Appendix 

El 
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Their expectation values are equal respectively to 

(*|^|*> = -ih , (4.4a) 

vdC e , </>=const 

(*||t)U|*) = Yl B ( v )\*( v , <l>)\ 2 • (4.4b) 

To calculate an explicit form of \& (needed to find the expectation values) one needs to 
find the spectrum of and the explicit form of its normalizable eigenf unctions. Because of 
the complicated structure of O, in order to do so one needs to resort to numerical methods. 
We present them in the next section. 

V. NUMERICAL STUDY 

This section is divided onto two parts. In section IV Al we present the methods and results 
of identifying the spectrum of the G operator and finding normalizable eigenfunctions. The 
techniques for computing the wave function and the expectation values are presented in 
section IV Bl In both parts, we applied the (appropriately refined) methods used already for 
the k = 1 model and introduced in [?J . Unless specified otherwise, from now on we will work 
with units in which G = 1. 

A. Spectrum of 

In order to construct the Hilbert space 7i phy e , one needs to find the eigenfunctions sup- 
ported on the lattice C e , which consists of two sublattices C±\ E \ '■= {±|e| + 4n; n e Z} 
invariant with respect to the action of the Hamiltonian constraint. Each of such eigenfunc- 
tions (denoted here as ip) is a solution to a difference equation: 

- u 2 B(v)ip(v) = C + (v)ij(v + 4) + (C°(v) + C A (v))ip(v) + C~(v)^(v - 4) , (5.1) 

where uo 2 is the eigenvalue that each given eigenfunction corresponds to and C°, C , C\ are 
given by (12.191 12.2ip . On each sublattice, this is a second-order equation - one needs to 
specify the initial data at two neighboring points (t>i n , v- m + 4) to uniquely define a solution. 
The symmetry condition ip(v) = tp(—v), however, restricts the amount of initial data in the 
following way: 

(i) For e G (0,2), the sublattices C±\ s \ are disjoint and the parity reflection II transforms 
one onto another. Therefore one needs to specify an initial data ip(v- m ), ip{v- in + 4) for 
just one of them, say £+\ e \, and complete it by the action of II. We denote such lattices 
as generic. 

(ii) When e = 0, 2 the sublattice £+\ e \ coincides with C~\ t \ and is invariant with respect 
to parity reflection. The condition ip(v) = ip(—v), applied to (15.11) . imposes on it an 
additional constraint of the form depending on the value of e: 

• 5 = 0: V(-4) = ip(0) = V>(4), 

• e = 2: if>(-2) = if>(2). Here the equality C~(2) = C + (-2) = implies addition- 
ally ^(±6) = -1(uj 2 B(2) + C A (2) + C°(2))/C + (2)]^(±2). 
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In consequence, the value of ip at just one point (v = or v = 2) determines the entire 
eigenfunction. These cases are called exceptional. 

The degrees of freedom specified above are complex; however, since the coefficients of 
( 15. ip are real, ip satisfies it iff so do its components $l(ip), ^(ip)- Therefore, we can safely 
restrict our study to a real ip. 

Upon this restriction, the space of solutions to (15. ip is 1-dimensional for exceptional 
lattices and 2-dimensional for generic ones. Once the initial data are specified appropriately 
for each case, the function ip can be found by solving (15. ip iteratively. 

To determine the properties of ip, we calculated the solutions in a wide range of both A 
([—10, — 1CT 6 ]) and uj ([0, 10 5 /i]). The qualitative features of the found solutions is visualized 
on figfSJ in general, for each ip one can distinguish 5 zones of distinct behavior, and the 
boundaries of these zones are specified by the functions vb(oj) and Vr(u), approximately 
equal to, respectively, the position of the bounce for a A = universe with = hu 
(determined in and the value of v at the recollapse point of the classical universe (given 
by (I3TT6D at <p = <jj" - <p ). 

(i) For \v\ < Vb(uj), the amplitude of ip grows/decays quasi-exponentially. 

(ii) For vb < \v\ < Vr, the behavior of ip is oscillatory (similar in nature to the behavior 

of dS3D). 

(iii) When \v\ > vr the eigenfunction grows/decays exponentially with |i>| (where the ex- 
ponential growth is a generic behavior). 

Note that for small uj, the zones (0} and (Jn]) may be empty (see figj3]for examples). 

Since we search for normalizable functions only, we have to select the ones which decay 
exponentially in the zones of type (lull) . We identify them numerically using different methods 
depending on whether the eigenfunctions are supported on generic or exceptional lattices. 

On exceptional lattices, each eigenfunction ip u is (for a given uj) determined uniquely up 
to a global scaling. To find the normalizable solutions, we scan the domain of uj using the 
following observation: 

Observation V.l For a chosen uj G [o^o^], ipu(e) — 1 an d v ^> vr(u> 2 ), the value ipuiy) 
is a continuous function of uj (more specifically, a polynomial) and its sign changes quasi- 
periodically . Furthermore, if we define to v ^ n as the values of u such that ^u, v>n (v) = 0, the 
limits uj n := lim^tx, io v%n are well defined and correspond to the values of u for which ij} u 
decays in zone (QUI) . 

In practice, due to the precision bound posed by numerical round-off, it is enough to 
(instead of finding the limits) look for values of uj v ^ n at v x sufficiently far away from Vr. For 
the actual search, we selected vt = max(2000, I.Svr). The search itself was performed in 
two steps: 

• First the sign of ^ w (vt) was checked for values of u> uniformly separated by a distance 
around 0.1. 

• If a change of sign was detected between neighboring points, the value of uj n ^ T was 
found via bisection. 
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FIG. 2: Examples of eigenfunctions of 6 supported on the lattices C £ for e = (a) and e = 1 (6). 

a) shows a normalizable eigenfunction of uj ~ 300.45 (red) and two divergent ones of a; respectively 
smaller (green) and larger (blue) by 0.1. For clarity, only the positive v part is shown. 

b) presents the absolute value of a normalizable eigenfunction of uj ~ 52.85 (red) along with two 
divergent examples: generic (green) and left-converging (blue) generated for, respectively uj ~ 53.35 
and 54.35. To show the behavior in a wide range of values, a logarythmic scale was used for the 
y-axis. 

Both figures correspond to A = —0.01. 



For generic lattices, the space of solutions is, up a to global rescaling, 1- dimensional, so 
besides uj we need to specify the value of ip at two points vi, v / + 4 e £+\e\- An additional 
complication is the fact that now the behavior in zones of type (lull) for v > and v < 
is independent. The function may grow for positive v while decaying for negative ones and 
vice versa. Therefore, to find the desired functions we divide the search procedure onto two 
steps: 



First we identify the family ifj w of functions decaying in zone (Imj) for v < (further 



denoted as left-converging) . To do so, we parametrize the initial data at Vj, v i + 4 by 
a parameter a G [0, ir] 

*l>aA v i) = cos(a) , ^a,w(«j + 4) = sin(a), (Q 

and scan the domain of a for the values at which the limit tim w _»_ oo , ct)W ('i>) = 0. 
Analogously to the exceptional lattice case, it is enough here to just choose some 
value — Vr(u)) 3> Vt- € £+[e| an d look for the values of a at which Tp a ,u( v T-) = 0. In 
practice, it suffices to choose u^- ~ —Vt, where vt is the value denned for exceptional 
lattices. The scan method is analogous to the scan of uj in the exceptional case: we 
divide the domain of a into 10 uniform intervals and if a change of sign of ipa tLU (vT-) 
is detected within an interval, the precise value of a is found via bisection. 

It was checked by inspection that, for each uj, there is exactly one value of a satisfying 
the above requirement. In consequence, for each uj the eigenspace of left-converging 
functions is 1-dimensional. 
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• Once the family ty w of left-converging functions is selected, we choose some vt ~ 
vt+ £ £+\e\ an d scan the domain of u for values at which ipu>{vT+) — 0, via the method 
specified for exceptional lattices. 




FIG. 3: The eigenfunctions eo to of 0, corresponding to A = —0.01 and supported on C £ with 
e = (a) and e = 1 (6). For clarity, only the v > part was shown in a) and only the part 
supported on C + \ E \ was shown in b). 



The search was first performed for small uo (< 50) to find the qualitative behavior of 
normalizable eigenfunctions. An example of the results is shown in Figs. [3] and Hh- All 
found eigenfunctions belong to one of the following groups: 

(1) Suppressed on the v < side with suppression exponential in uj. 

(2) Suppressed for v > 0. 

(3) Peaked about v — 0. 

In consequence, it is most convenient, from the point of view of the numerical precision 
of the solutions, to specify the initial data at vj ~ ±^/j. However, because of the quasi- 
exponential behavior of the eigenfunctions in zone (OJ, we can calculate (with a sufficiently 
small numerical error) only the solutions suppressed on the side where the initial data were 
specified. Therefore it is necessary to repeat the search twice: for vj ~ vr and vj ~ — Vr. 

The spectrum scan described above was performed for 18 values of A ranging from —20 
to — 10~ 6 . It revealed the following properties (visualized in figs. [3] -[5]). 

• As analytically predicted, for each e the spectrum of O is discrete and the eigenvalues 
are isolated. With the exception of the lowest u, the eigenfunctions are highly (ex- 
ponentially in y/uj) suppressed for one triad orientation (sign of v). The eigenvalues 
corresponding to them are continuous functions of e. The density of eigenvalues is 
twice higher on the generic lattices than on the exceptional ones. Furthermore, for 
e = 2, the two families of left-suppressed and right-suppressed eigenfunctions converge 
(see fig. Sh). 
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FIG. 4: a) The lowest (cj < 44) elements of 0's spectrum are shown as functions of ±|e|. The 
eigenvalues are divided into three groups corresponding to the following eigenfuctions: (1) left- 
suppressed (red crosses), for which > 50% of the norm is located on v > 0, (2) right-suppressed 
(green x-es) defined analogously and (3) singularity-peaked (blue stars), where > 50% of the norm 
is located at the three points closest to v = 0. 

b) The large u limit of the eigenvalue separation Alo, shown as function of A (red crosses). The 
blue line represents the small A limit given by (|5.5p . 



• The separation Au n := u> n+ i — u n is not uniform. It depends on e and A as well as n. 
However, for large values of a;, Au n converges to the limit value Auj with convergence 
rate u~ 2 (see ngJSH) 

Au n = Au + 0{u~ 2 ) , Auj = lim Au n . (5.3) 

n^oo 

Numerical inspection shows that the correction satisfies (with the exception of the 
lowest a;) the following bound relation 

\Auj-Auj n \ < (5.4) 
or 

where, for |A| < 10, A < 0.21 and A decreases for smaller |A|, reaching in the |A| — ► 
limit the value A w 0.1358 ± 2 ■ 10" 4 (see figEJ). 

• The limit Ao; was found numerically via 4th order polynomial extrapolation of Au n . 
It is a function of A only, i.e. it does not depend on e. Its values for different superse- 
lection sectors agree up to 10~ 9 precision. The dependence on A found numerically is 
shown on fig. H)d. For small values of |A| it can be approximated via a power function 

Au^a\AG\ b (5.5) 

where a ta 3.87 and b w 0.0489. 

The spectrum and normalizable eigenfunctions found here may be next used to construct 
the semiclassical states. Details of this construction will be presented in the next section. 
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FIG. 5: The rescaled eigenvalue separation correction term A(u) := \Au n — Au\(ui n / Au>) 2 is shown 
in a) for several values of A. Its uj — > oo limit is plotted in b) as a function of A. The 'wiggles' at 
small values of | A| are the results of numerical errors due to the precision limitation of the applied 
calculation method. 



B. Semiclassical states, evolution 

Once we know the values of u) n and e n (v), the construction of a physical state from (14.11) 
is straightforward. There are two possibilities here: direct summation of equation (14.1b ) or 
numerical integration via equation (14.1b ) (or equivalently via ( 12.181) ) of some initial data 
specified at a given <f) . To find these data, we again have two methods at our disposal: 
one of them is the same direct summation of (14.1b ). but applied to one slice, whereas the 
second possibility is the use of a slice of a WDW semiclassical state (see section HTT|) peaked 
at large v*, where we do not expect strong quantum-geometric effects. In practice we used 
the second method, integrating the state in via equation (12 . 1 8[) and using as initial data 
both the WDW slices and the results of the summation of (14.1b ). The first method of state 
calculation was used only to measure the wave packet spread increase in large intervals of 
<ft, as the integration methods were not precise enough for this application. 



1. Initial data 

Let us focus on the second method of initial data specification: constructing the WDW 
slice. In order to be able to directly compare the dynamics of LQC model and its WDW 
limit described in section II III we take as the initial data the = O section of the Gaussian 
state (13.141) peaked at = ftw* and v*. Since (12.181) is a second order equation, to specify 

the initial data completely we also need \l/(f , (f) ) - the first order derivative of \1/ with respect 
to 4>. We get it by integrating the integrand of (I3.14p . multiplied by iui, over to. 

In order to calculate the specified integrals, we first need to compute the values of e^iy), 
which are the normalized Bessel functions K, (see section UlIBj) . To do so, we apply the 
combined method specified by Gil, Segura and Temme [201 ]. 

Once we have e^(f), we integrate ( 13.141) (and the analogous expression for over the 
domain [u* — 7 a, u* + 7 a], using the trapezoid method. Such choice of domain provides 
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sufficient precision - the errors due to the removed tails are much smaller than the error 
associated with the computation of the K, functions. 

Note that we intend to construct the initial data corresponding to the positive frequency 
solution to (I2.18p . In that case, * is already determined by \l/ via (14.1a ). On the other 
hand, we determined it using positive frequency WDW equation (13.91) . Since y/Q differs 
from our initial data is not a pure positive frequency solution. To minimize the negative 
frequency part, we choose the following method to construct states sharply peaked at large 
v*: we require v* be greater than 2.5p^/h, which keeps the negative frequency part below 
10 -3 of the entire wave packet norm. 

We avoid the above problem if we use directly the basis of functions e n (v ) and sum them 
using (14.1b ). In that the spectral profile \l/ n we choose the restriction of the Gaussian 

to {u n }, that is 

m n = e sa— e -1 ^* (5.6) 

where ftiu* = is again the location of the peak in the momentum and fkr/y/2 is its spread. 
The parameter 0* is determined by the position v* of the peak in v and value of <p via 
(I3.15p . Similarly to the WDW initial data, we sum only over u n G [oo* — 7cr, uj* + 7a]. The 
derivative * is calculated by summing over the individual terms, multiplied by iu n . 



2. Evolution 

Given some initial data, one can integrate it over some interval [<f) , 0i] using equation 
(I2.18p . which is a system of a countable number of coupled ordinary differential equations 
(ODE). Due to the v -reflection symmetry, it is enough to restrict the domain of integration 
to C + = C+\ v \ for generic lattices and C + = C e PI 1R + for exceptional ones. Additional, 
the numerical nature of our study requires that we further restrict the domain of v to the 
finite subset £^ nax := C + PI [— v max , f max ], imposing at the outermost points of the domain 
some (artificial) boundary conditions. Since the system under consideration is a classically 
recollapsing one, it is enough to choose the reflective conditions ^ = ^ = 0. To prevent their 
interference with the dynamics, we have chosen t> max to be not smaller than 1.3vji(u!*) + 2000. 

Upon the above restriction of the v domain, the equation (12.181) becomes a finite system 
of ODEs. We integrate it using a Ath-order adaptive Runge-Kutta method (RK4). To adapt 
the steps of integration, we compare solutions corresponding to step A0 and A0/2 and 
require the difference between them (at a single A0 step / two A0/2 steps) to satisfy the 
inequality: 

- *A0/2 < H ^11^/211, (5.7) 

\<Pl-<Po\ 

where e is a preset global bound. The two solutions are compared via the following norm 

H^ll := sup \ty(v,(f))\ . (5.8) 

Since only |^| enters the formulae for the expectation values of \v\,p and u|, it is also conve- 
nient to introduce an auxiliary metric measuring the error in absolute values only 

||*i -* 2 |U = SU P | l*i(^ </>)!- I*2(^,0)|| (5.9) 
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An example of convergence test done with respect to both the norm (15.81) and the metric 
( 15. 9p is shown in fig® where the results of integration with different error bounds e were 
compared against the result of polynomial extrapolation at e = 0. As we can see, the 
integration error is for at least one order of magnitude smaller than that for \1/ itself. 




10 4 10 5 10 6 10 7 

N 



FIG. 6: Convergence test for the integration method of a Gaussian wave packet generated with 
A = —0.1 and peaked at p^ = 10 3 , with relative p^ spread 0.05 and v* = 0.5vr(p^). The initial 
data were specified at 4> = and evolved till = 1. The upper (red) curve shows the norm 
of the difference H^at) — ^|| between the slice (j) = 1 of the solution ^(tv) corresponding to the 
integration with N steps and the same slice of its N — > oo limit $ (found via 8th order polynomial 
extrapolation). The lower (green) curve shows the analogous difference taken with respect to the 
metric (|5.9p . 



3. Observables 

Knowing an explicit form of $ at £^~ max x [<f) , cpi], we can complete it to (£^ max U C~ ) x 
[4>oi <fii) (where £~ max := {—v : v G }) via reflection and find the expectation values of 
the observables (14. 31) . restricting the sums (14.21 |4.4p to a finite domain £+ U C~ . Their 
dispersions can be in turn calculated in the standard way 

<ak>u> 2 = (vi)-(\v\,) 2 (A^f = (%) -(^f . dsnn 

where ( ) , ( (fy ) have a form analogous to (I4.4p . 

In addition to \v\^, (p^, it is useful to introduce another family of observables: the regu- 
larized energy density at a given moment of </> 

(-M ptB+ptBt, := e^-^B(v)^>(v,<P') . (J5H) 



We calculate their expectation values via 



w*i*> = S b{v)\<$>{va)\\ $ = w^*), (5.i2) 

pi V 7/ c+v uc - v 



m ax <-* *— %ax 
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whereas the dispersions can be derived analogously to (I5.10p . 

The above methods for calculating the expectation values were applied to the wave func- 
tions calculated earlier through the RK4 method. We analyzed the states evolved (inte- 
grated) from both WDW and exact LQC Gaussian wave packets corresponding to 17 values 
of A ranging from —20 to — 10~ 6 , for 5 different superselection sectors covering the full range 
of e. The peak in momentum 0^ covered the values from 5 • 10 2 to 10 4 (10 values), while its 
relative spread ranged from 0.01 to 0.1. 

VI. RESULTS AND DISCUSSION 




FIG. 7: The absolute value of the wave function representing a Gaussian state (|4.ip generated via 
backward integration of an initial profile corresponding to A = —1, = 5 • 10 3 H, Ap^/p^ = 0.01, 
v* = 0.6 vjiip^) and evaluated at 4> Q = 0. For presentation clarity, only values > 10~ 6 were shown 
on the plot. 

An example of the results of our numerical investigations is presented in figs. [7]-[TUJ The 
general properties of the considered model are similar to the ones of the models previously 
investigated: A = and k = l, that is: 

• The states remain sharply peaked for long evolution times. On each superselection 
sector and large u, the spectrum of G quickly approaches uniformity (with approach 
rate u~ 2 ). In consequence, a wave packet sharply peaked at large should be almost 
periodic in (f). This expectation is confirmed by our numerical results, where already 
for of the order of few thousands the departures from periodicity were undetectable 
within given precision of integration. 
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• For large volumes (small energy densities), the trajectory of the expectation values 
(H</>) agrees with the classical one given by (I3.16p . In particular, the universe recol- 
lapses at the volume predicted by the classical theory even for large values of A; this 
was numerically confirmed for |A| up to 20. 

• Once the expectation value of the energy density approaches the Planck order, we ob- 
serve the departures from the classical theory due to quantum-geometric corrections. 
The corrections act effectively like an additional repulsive force, which in particu- 
lar causes the bounce at the point where the total energy density (p<f,) + A/(87rG) 
approaches a critical value p c ~ 0.82ppi, identified already in 

• After the bounce, the universe again enters (another) classical trajectory repeating 
the cycle of expansion, recollapse and contraction till the energy density grows again 
to Planck scale. In consequence, the evolution is periodic and, similarly to the k = 1 
case, we are dealing with a cyclic model. 

• The wave packet remains sharply peaked even in the region where the quantum cor- 
rections are strong. In consequence, the evolution can be described by the classical 
effective dynamics, similarly to the A = case. Indeed, the comparison of the values 
of ( \v\cj)) with the effective trajectories corresponding to the holonomy corrections (see 
Appendix IB ip has shown that they agree up to an error well below (Av). 
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FIG. 8: The expectation values (red bars) of (o) and p^ (b) are compared against the clas- 
sical (red lines) and effective (blue lines) trajectories of v((ft) and p<j>{v) respectively. The data 
corresponds to a Gaussian wave packet (|4.ip with A = —0.1, p^ = 10 K, Ap^/p^ = 0.012, 
v* = 0.55 vji(p^ ) ) specified at cj) = and evolved backwards. Because of the large changes in 
magnitude of (p^), a logarythmic scale was used for the y-axis of fig. b). 



The results listed above show that the picture based on the analysis of the previous 
models is valid here as well. Similarly to that cases, the correction due to the discreteness of 
geometry cause gravity to become repulsive at large energy densities and, in particular, force 
the universe to bounce when the energy density reaches a critical value. This indicates that p c 
may be a fundamental quantity, independent on the matter content at least in isotropic cases. 
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FIG. 9: A detailed picture of the comparizon of (H</,) against the v(4>) trajectories presented in 
fig. [Bh, is shown near the bounce (a) and recollapse (6) points respectively. 



1 

0.9 
0.8 
0.7 
0.6 
0.5 
0.4 
0.3 
0.2 



a) 









LQC 
classical 
effective 




1 

i 


1 

i 

I 




,1 

/ f 


1 

I 

I 


s. 

E 

3 

i 


\ \ 


/ f 
/ I 

X 

/ 1 

¥ 

1 

i 

t' 






\ \ 
I. \ 
'i. \ 
1. 

1. 

?, 

S. ' - 

s. ■ 
s.\ 
% 

3 



0.005 
P 

0.0045 



b) 



LQC 
classical 
effective 



-3.3 -3.25 -3.2 -3.15 -3.1 -3.05 -3 -2.95 -2.9 



-2.55 -2.5 -2.45 -2.4 -2.35 -2.3 



FIG. 10: A detailed picture of the comparizon of (p^) against the p^cj)) trajectories presented in 
fig. [8b is shown near the bounce (a) and recollapse (b) points respectively. 



Furthermore, the states remain sharply peaked even in regions where quantum-geometric 
effects dominate the dynamics, where in principle one expects to loose the semiclassicality. 
The dynamics itself can be well approximated by an effective Friedmann equation (see 
Appendix IB ip 

H 2 = «fp(l - p/p c ) , (6.1) 

where if is a Hubble rate and p is a total energy density. 

The agreement between the quantum evolution and the effective one brings out another 
issue: since the spectrum of is not exactly uniform, the states are not exactly periodic 
and a spread increase can be observed between cycles. This leads ultimately to the loss of 
semiclassicality. This in turn raises the question about the size of time interval in which the 
state remains sharply peaked. 

To answer this question, we analyzed the spread increase within one cycle of evolution. 
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It can be estimated via the heuristic methods described in section IB 21 and turns out to be 
much smaller than in k = 1 case. For example, when A m — 1CT 120 , a universe peaked about 
p~0 large enough for it to grow to megaparsec size, and with relative dispersions in p^ and v 
of the same order, will need at least 10 70 cycles for the relative dispersion to double. The 
number of cycles needed to grow to a considerably large value (say 10~ 6 ) is correspondingly 
larger. 

For small values of the momentum (that is pi < 10 3 h), we were able to confirm the 
heuristics numerically. Also, since for larger momenta the states become more and more 
semiclassical, we expect the estimate to become more accurate there. The result is however 
far from general, as numerical tests were done for a specific family of states only, thus 



(as it was discussed in [211 ]) do not allow us to exclude the situation, where some specific 
example of state violates the bound. On the other hand, the proposed estimate is based 
on the properties of the spectrum of 9, thus we expect that a bound of at least a similar 
order should hold in general. Such situation happened for example in the A = case fl~lj ]. 
where it was possible to find (in the context of sLQC) an analogous bound satisfied by all 
the states which admit semiclassical epoch (see section [T|) in their history. A similar bound 
was next derived in exact LQC [22|. More precise statements regarding model considered 
here will however require further work. 

Acknowledgments We would like to thank Abhay Ashtekar, Wojciech Kamihski and 
Parampreet Singh for extensive discussions and helpful comments. We also profitted from 
discussions with Martin Bojowald and Jerzy Lewandowski. This work was supported in 
part by the National Science Foundation (NSF) grant PHY-0456913 and the Eberly re- 
search funds. TP acknowledges financial aid provided by the I3P framework of CSIC and 
the European Social Fund. EB acknowledges the support of the Center for Gravitational 
Wave Physics, funded by the National Science Foundation under Cooperative Agreement 
PHY-01 14375. 



APPENDIX A: ANTI-SYMMETRIC SECTOR OF LQC 

Due to the lack of a parity violating interaction in the model considered in this article, 
the change in the triad orientation is a large gauge symmetry. This allowed us to restrict the 
physical Hilbert space to the subspace of states invariant with respect to the reflection in v 
corresponding to this orientation change - the symmetric sector. In principle, however, we 
could choose instead the space of states which are antisymmetric under considered transfor- 
mation. There is no physical reason to favor one of these two choices over the other. This 
raised a concerns on whether the results of LQC are tied to the selection of symmetric sector 
and whether they will still hold in the antisymmetric one. We address these concerns here 
by repeating the constructions of section [TVl this time building the physical Hilbert space 
out of antisymmetric states. 

First, following section [TV] we divide the kinematical Hilbert space 7igrav onto supers- 
election sectors, i.e. the spaces 'H^ v t of functions supported on lattices C £ . The results 
of 0] (self-adjointness of and discreteness of its spectrum on each of these spaces) were 
derived without any symmetry assumption, thus they hold also in our case. Furthermore, 
as we will show below, the spectrum is non-degenerate also when we restrict the space of 
eigenfunctions to the antisymmetric ones. In consequence, we can construct the physical 
Hilbert space as specified in (14.11) . but by imposing on the relevant eigenfunctions the 
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condition e"(t>) = — e^(— v) instead of the symmetry requirement. 

To check the effect of the above modification on the dynamics, we have to examine how 
it changes the exact form of e n . That, in turn, depends on the value of the superselection 
sector label e. 

• For e 7^ 0, 2 (generic lattices), the symmetric eigenfunction on C e is completely deter- 
mined (see discussion in section IV Ap by its restriction to the lattice C + \ e \, with the 
remaining part supported on £_i e i determined via a symmetry relation. Furthermore, 
symmetry does not impose any constraint on the part supported on C+\ e \ itself and we 
can complete it to the antisymmetric eigenfunction by simply acting with —II on it. 
In consequence, there exists a 1 — 1 correspondence between these two types of eigen- 
functions. Namely, each antisymmetric eigenfunction ip a is related to the symmetric 
one if> y i a: 



This implies that, in the antisymmetric sector, the spectrum of G is the same as in 
the symmetric one. 

• When e = 2, the situation is similar to the generic case. The solutions to (15.11) on two 
sublattices C £ fl M + and C £ PI K~ are independent, thus each eigenfunction is again 
determined by its restriction to C £ fl K + . In consequence, we again have the 1 — 1 
correspondence between the symmetric and antisymmetric eigenfunctions 



and the spectra of in both sectors are identical. 

• The case when e = requires a bit more care. In section IV A| the symmetry assumption 
imposed on solutions to (15.11) an additional constraint, allowing to determine if} (4) for 
known if>(0). Antisymmetry replaces this constraint by a different one: ip a (0) = 0. 
Therefore the whole procedure of identifying the normalizable eigenfunctions has to 
be redone. We can however apply exactly the same procedure as in section IV Al The 
results are as follows: 

— The qualitative features of the eigenfunctions remain the same. In particular, 
we can still distinguish the same 5 zones of exponential/oscillatory behavior (see 
figUTa). Their boundaries vb, Vr are exactly the same as in the symmetric case. 

— The spectrum of in the antisymmetric sector is different than in the symmetric 
case, however the eigenvalues of one sector approach the ones of the other very 
quickly (see figJTTb). In consequence, the separation between the eigenvalues 
approaches, as u — > oo, the same limit shown in fig. Hb. The rate of approach to 
this limit also remains the same. 

The similarity between eigenfunctions of the two considered sectors implies an analogous 
similarity between the physical states. In particular, for e ^ 0, if the eigenfunctions satisfy 
equations (IA1HA2|) . so will the wave functions. Then if we take two physical states, a 
symmetric and an antisymmetric one with the same spectral profile ty n , the expectation 




(Al) 




(A2) 
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values of (all the powers of) the observables defined in section IV B 31 will be exactly the same 
for both of them. 

For e = 0, due to the slight difference in the spectrum, we have to repeat the analysis 
of IV Bl But again the numerical checks reveal no measurable deviations from the results 
obtained in the symmetric sector. 

In summary the results obtained for both e 7^ and e = show that working with the 
antisymmetric sector instead of the symmetric one does not produce any qualitative changes 
or (apart from a slightly different spectrum of in e = case) any measurable modifications 
to the physics predicted by the model. 
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FIG. 11: a) A set of lowest (a> < 44) spectrum elements of for the symmetric and antisymmetric 
sector (with A = —1) is shown with respect to ±|e|. The green x-es represent the eigenvalues cor- 
responding to the cases where the relation between symmetric and antisymmetric eigenfunctions 
is given by (|A11 IA2[) (denoted as generic) . The red crosses and blue stars represent the eigenvalues 
of, respectively, the symmetric and antisymmetric eigenfunctions on the lattice £ e =o- The anti- 
symmetric eigenfunctions eg to e\ corresponding to the eigenvalues shown in a) are presented in 
b). For clarity, they are plotted on the v > semiaxis only. 



APPENDIX B: HEURISTIC DESCRIPTION 

In this sections we discuss some issues related to the heuristic method for the description 
of the quantum evolution. We divide its content into two parts, dedicated respectvely to the 
effective classical dynamics and the estimate of the dispersion growth during the evolution 
of the semiclassical state. 



1. Effective dynamics 

The numerical tests described in the main body of the paper have shown that if a state is 
semiclassical at some epoch, it will remains so for a large fraction of the evolution (i.e., a large 
number of cycles of bounces and recollapses) . In particular it remains sharply peaked even 
in the regions where the quantum gravity corrections modify the dynamics. This indicates 
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the existence of a (n effective) classical theory whose predictions well agree with those of 
LQC. 

At the rigorous level, such theory was derived for A = [23| with the use of the geo- 
metric formulation of quantum mechanics 24j. Up to the second order quantum corrections 
(remaining always small during the evolution}, its results confirm the predictions of the 
classical effective dynamics proposed earlier [HI, [6|, derived heuristically by replacing the 
connection c in classical Hamiltonian by sin(/ic)//i. Here we apply this heuristic method 
to the system with a cosmological constant considered in the main body of the paper. An 
analogous derivation of the effective dynamics (to the level of quadratures) and the analysis 
of the resulting trajectories was done in [l5| . however the trajectory parametrization used 
there makes the direct comparizon with the results of quantum evolution difficult. 

The classical Hamiltonian of the system is related to the constraint (12. 4[ 12.61) via 7i c g = 
C/(16nG). Application of the rule c — > sin(/ic)//2 yields the result 

Q 2 - 

8n&Yp 2 p j 167tG 

Hamilton's equations v = {v, H e s} and <fi = {(ft, Ti. e s} are identical to the A = case. Written 
in terms of v they are respectively 

Sv , ( 6 \% K 

sm(yuc) cos(yuc) 
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Taking the square of (1B2a ) and supplying sin(/2c) via ( 1B1I) . we arrive to an analog of Fried- 
mann equation: 

where p and p c are the total matter energy density and the critical energy density found in 

a 

A p 2 a/3 

Applying flB2b ). we can rewrite the resulting Friedmann equation in the form involving v 
and <p only 



Vcj> 



± W12ttG 



(B5) 



The sign in front of the right hand side of the above equation depends on the evolution 
epoch, and in particular changes during recollapse. Therefore, it is convenient to rewrite 
(1B5[) in the second order form (obtained by differentiating it): 

v H = 12nGv \(^ -l)(l-P-) + P-} . (B6) 

AP* J \ Pc) Pc. 

To compare the results of the numerical evolution, we integrated equation (1B6I) numeri- 
cally using a fifth-order adaptive Runge-Kutta method (known as RKJ^S). The initial value 
v, needed to complete the initial data specification, was calculated via (1B5I) . 

An example of the comparison results is shown in figs. UHl The trajectories agree 
with the results of the quantum evolution everywhere. The differences between them are 
much smaller than the spreads of the wave packets even near the bounce. 
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2. Bound on the dispersion growth 



The analysis of section |V] has shown that the states that are semiclassical at a given 
initial time O remain so for many cycles of bounces and recollapses. However, due to 
non-uniformity of the spectrum of the operator, the initially coherent wave packet slowly 
spreads out. Here we derive an upper bound on this spread growth using some heuristic 
estimates and applying the knowledge about the spectrum of presented in section IV Al 

To start with, let us note that for large p^ = fruj the distance between neighboring 
eigenvalues is almost constant uj n +i — u n ~ Au (see ( 15. 3p ). In consequence, the wave 
function is almost periodic in 0, with an approximate period equal to 

T "^- < B7 > 

Now, if we consider two classical (effective) trajectories corresponding to p^ equal respec- 
tively to fojj and ftw + 5u, the difference between periods is determined by the corrections to 
the uniformity of Au n . They are in turn bounded by the function A (Alj) 2 oj~ 2 (see ( 15.41) ). 
Applying this bound to (1B7|) (i.e. taking Au n = Au(l + Auj~ 2 )) and neglecting terms of 
higher order in 8u, we obtain the following difference in T: 

4-7T A 

5T w —5u , (B8) 

to 6 

which can be now used to estimate the growth of Av/v within one cycle. To do so, we note 
that, since the cosmological constant term acts like a positive v 2 potential, the speed is 
bounded from above by the speed of a classical universe with A = 

|^| < \v;\ := V12ttG\v\ . (B9) 

In consequence: 

S ^<8^A\^. (BIO) 

In order to arrive to this bound, we used some heuristic methods that need to be confirmed 
using numerics. Unfortunately, due to the extremely small values of 8v/v, we were able to 
check (IBlOp only for small values of the frequency, u> < 10 3 . To do so, we calculated the 
semiclassical states in two intervals of (ft separated by a large (> 100) distance in (ft and 
compared the difference between the maximal relative dispersion in v observed within one 
cycle in both of the chosen intervals. To compute the wave functions, we used a direct 
summation method specified in section IV B[ Within the checked range 500 < uj < 1000, the 
bound was satisfied. 
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